# Define constant 10x smaller than minimum nonzero biomass
const <- min(data_kelp$kg_per_m2[data_kelp$kg_per_m2 > 0])/10
data_sp <- data_kelp %>%
mutate(
# Add small constant
kg_per_m2_c = kg_per_m2 + const,
# Log-transform to variable with small constant
log_kg_per_m2_c = log(kg_per_m2_c))
# Select term columns to join with residuals df
data_terms <- data_sp %>%
dplyr::select(obs_id, age_at_survey, site_type,
all_of(grep("^(hard|soft|kelp|depth)", names(.), value = TRUE)))
data_sp_scaled <- data_sp %>%
mutate(across(where(is.numeric), ~ scale(.) %>% as.vector()))Model Diagnostic Comparisons
Overview
Diagnostics from initial models demonstrated moderate issues with kurtosis, nonlinearity, and heteroskedasticity. The goal of this comparison is to explore whether other modeling options improve these issues. The first focal species is vermillion rockfish (Sebastes miniatus), a species found statewide that is generally associated with hard bottom habitat. This was chosen because it has two int
The original approach was:
Generalized linear mixed effects models with all numeric variables scaled (including dep. variable, first log+1 transformed then scaled to mean = 0 and sd = 1). The variables were scaled to accommodate differences in scales between MPA age and the habitat variables, and the magnitude between the various buffer scales now that multiple scales could be included in the same models (e.g. hard bottom at 500m and kelp at 25m, varying from 10s of sq meters to millions of sq meters).
The top model for this species was:
Log+1-Transformed Biomass ~ Hard Bottom 500m * Site Type + Kelp Annual 100m * Site Type + Site Type * MPA Age
- Does log-transforming biomass make a difference in model fit? Is there a difference between using no transformation, log+1, vs. log+c, where c is a very small constant (10x smaller than the smallest nonzero vale)
- Does scaling the variables (mean = 0 and sd = 1) make a difference in conclusions or model fit?
- Does adding a 3-way interaction with MPA age make a difference in model fit?Build
Untransformed
B ~ H*ST + K*ST + ST*A (Not Scaled)
- Removed year as it had zero variance (singular fit)
- Significant:
H*ST, K*ST, ST*A, ST
glm1 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescalingattr(glm1, "name") <- "GLMM: B ~ H*ST + K*ST + ST*A (Not Scaled)"
summary(glm1)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula:
kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type +
site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
Data: data_sp
AIC BIC logLik deviance df.resid
-15536.1 -15475.9 7779.0 -15558.1 1737
Scaled residuals:
Min 1Q Median 3Q Max
-2.6605 -0.3363 -0.0311 0.1542 14.3967
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 9.310e-07 0.0009649
bioregion (Intercept) 1.968e-06 0.0014028
Residual 7.744e-06 0.0027829
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.975e-03 9.058e-04 3.631e+00 2.180 0.101628
hard_bottom_500 -9.557e-10 5.811e-10 1.630e+03 -1.645 0.100214
site_typeMPA -3.139e-03 4.155e-04 1.721e+03 -7.555 6.75e-14
kelp_annual_100 -3.892e-08 2.924e-08 1.745e+03 -1.331 0.183258
age_at_survey 1.928e-05 1.785e-05 1.743e+03 1.080 0.280182
hard_bottom_500:site_typeMPA 1.009e-08 1.080e-09 1.572e+03 9.342 < 2e-16
site_typeMPA:kelp_annual_100 1.529e-07 3.639e-08 1.745e+03 4.203 2.77e-05
site_typeMPA:age_at_survey 9.706e-05 2.567e-05 1.734e+03 3.781 0.000161
(Intercept)
hard_bottom_500
site_typeMPA ***
kelp_annual_100
age_at_survey
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.234
site_typMPA -0.146 0.297
klp_nnl_100 -0.124 0.130 0.197
age_at_srvy -0.185 0.061 0.361 0.246
h__500:_MPA 0.041 -0.309 -0.740 -0.003 0.033
s_MPA:__100 0.099 -0.124 -0.255 -0.732 -0.200 -0.037
st_tyMPA:__ 0.139 -0.074 -0.521 -0.161 -0.681 -0.050 0.215
fit warnings:
Some predictor variables are on very different scales: consider rescalingB ~ H*ST + K*ST + ST*A (Scaled)
- Removed year as it had zero variance (singular fit)
- Significant:
H*ST, K*ST, ST*A, ST - Same as nonscaled version (#1)
glm2 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)
attr(glm2, "name") <- "GLMM: B ~ H*ST + K*ST + ST*A (Scaled)"
summary(glm2)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula:
kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type +
site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
Data: data_sp_scaled
AIC BIC logLik deviance df.resid
4306.8 4366.9 -2142.4 4284.8 1737
Scaled residuals:
Min 1Q Median 3Q Max
-2.6605 -0.3363 -0.0311 0.1542 14.3967
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 0.07924 0.2815
bioregion (Intercept) 0.16748 0.4092
Residual 0.65916 0.8119
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.19558 0.25263 3.04429 0.774 0.494481
hard_bottom_500 -0.04452 0.02707 1630.32168 -1.645 0.100214
site_typeMPA 0.28866 0.04080 1742.29575 7.075 2.16e-12
kelp_annual_100 -0.04514 0.03390 1744.51016 -1.331 0.183258
age_at_survey 0.03029 0.02804 1743.12071 1.080 0.280182
hard_bottom_500:site_typeMPA 0.47000 0.05031 1571.97148 9.342 < 2e-16
site_typeMPA:kelp_annual_100 0.17734 0.04220 1745.12473 4.203 2.77e-05
site_typeMPA:age_at_survey 0.15247 0.04032 1734.11137 3.781 0.000161
(Intercept)
hard_bottom_500
site_typeMPA ***
kelp_annual_100
age_at_survey
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.024
site_typMPA -0.080 -0.007
klp_nnl_100 0.003 0.130 -0.030
age_at_srvy 0.007 0.061 0.008 0.246
h__500:_MPA -0.013 -0.309 0.026 -0.003 0.033
s_MPA:__100 -0.003 -0.124 -0.021 -0.732 -0.200 -0.037
st_tyMPA:__ 0.003 -0.074 -0.024 -0.161 -0.681 -0.050 0.215 B ~ H*ST*A + K*ST*A (Not Scaled)
- Removed year as it had zero variance (singular fit)
- Significant:
H*ST*A, K*ST, ST*A - Almost significant:
H, ST
glm3 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescalingattr(glm3, "name") <- "GLMM: B ~ H*ST*A + K*ST*A (Not Scaled)"
summary(glm3)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula:
kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 *
site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
Data: data_sp
AIC BIC logLik deviance df.resid
-15611.3 -15529.3 7820.6 -15641.3 1733
Scaled residuals:
Min 1Q Median 3Q Max
-2.9976 -0.3236 -0.0405 0.1645 14.7939
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 9.166e-07 0.0009574
bioregion (Intercept) 2.114e-06 0.0014539
Residual 7.381e-06 0.0027167
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df
(Intercept) 2.155e-03 9.572e-04 3.952e+00
hard_bottom_500 -1.674e-09 9.524e-10 1.744e+03
site_typeMPA 4.896e-04 6.226e-04 1.739e+03
age_at_survey -2.231e-06 3.181e-05 1.740e+03
kelp_annual_100 -2.060e-08 4.054e-08 1.746e+03
hard_bottom_500:site_typeMPA -1.416e-09 1.875e-09 1.717e+03
hard_bottom_500:age_at_survey 8.609e-11 9.125e-11 1.736e+03
site_typeMPA:age_at_survey -3.415e-04 6.185e-05 1.737e+03
site_typeMPA:kelp_annual_100 1.342e-07 5.380e-08 1.737e+03
age_at_survey:kelp_annual_100 -3.951e-09 5.866e-09 1.744e+03
hard_bottom_500:site_typeMPA:age_at_survey 1.380e-09 1.870e-10 1.739e+03
site_typeMPA:age_at_survey:kelp_annual_100 2.683e-09 7.479e-09 1.731e+03
t value Pr(>|t|)
(Intercept) 2.251 0.0884 .
hard_bottom_500 -1.758 0.0790 .
site_typeMPA 0.786 0.4318
age_at_survey -0.070 0.9441
kelp_annual_100 -0.508 0.6114
hard_bottom_500:site_typeMPA -0.755 0.4504
hard_bottom_500:age_at_survey 0.943 0.3456
site_typeMPA:age_at_survey -5.522 3.86e-08 ***
site_typeMPA:kelp_annual_100 2.494 0.0127 *
age_at_survey:kelp_annual_100 -0.674 0.5007
hard_bottom_500:site_typeMPA:age_at_survey 7.379 2.45e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100 0.359 0.7199
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500 -0.316
site_typMPA -0.175 0.398
age_at_srvy -0.292 0.669 0.429
klp_nnl_100 -0.130 0.118 0.185 0.309
hr__500:_MPA 0.119 -0.431 -0.889 -0.318 -0.053
hrd__500:__ 0.234 -0.803 -0.353 -0.813 -0.094 0.405
st_tyMPA:__ 0.154 -0.350 -0.830 -0.501 -0.155 0.736 0.413
s_MPA:__100 0.094 -0.088 -0.179 -0.229 -0.714 -0.040 0.068
ag_t_:__100 0.054 -0.035 -0.108 -0.248 -0.708 0.043 0.065
h__500:_MPA: -0.116 0.391 0.740 0.387 0.050 -0.826 -0.481
s_MPA:__:__ -0.037 0.018 0.069 0.187 0.536 0.058 -0.050
st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500
site_typMPA
age_at_srvy
klp_nnl_100
hr__500:_MPA
hrd__500:__
st_tyMPA:__
s_MPA:__100 0.122
ag_t_:__100 0.123 0.530
h__500:_MPA: -0.897 0.065 -0.033
s_MPA:__:__ -0.060 -0.748 -0.761 -0.113
fit warnings:
Some predictor variables are on very different scales: consider rescalingB ~ H*ST*A + K*ST*A (Scaled)
- Removed year as it had zero variance (singular fit)
- Significant:
H*ST*A, H*ST, K*ST, ST*A, ST - Almost significant:
H - Similar to non-scaled version except site type is now highly significant on its own
glm4 <- lmer(kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)
attr(glm4, "name") <- "GLMM: B ~ H*ST*A + K*ST*A (Scaled)"
summary(glm4)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula:
kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 *
site_type * age_at_survey + (1 | bioregion) + (1 | affiliated_mpa)
Data: data_sp_scaled
AIC BIC logLik deviance df.resid
4231.6 4313.6 -2100.8 4201.6 1733
Scaled residuals:
Min 1Q Median 3Q Max
-2.9976 -0.3236 -0.0405 0.1645 14.7939
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 0.07802 0.2793
bioregion (Intercept) 0.17991 0.4242
Residual 0.62820 0.7926
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df
(Intercept) 0.18751 0.26058 3.00887
hard_bottom_500 -0.04453 0.02646 1635.65954
site_typeMPA 0.28244 0.04112 1741.50636
age_at_survey 0.02530 0.02866 1743.68488
kelp_annual_100 -0.06210 0.04063 1743.73621
hard_bottom_500:site_typeMPA 0.47020 0.04928 1585.04419
hard_bottom_500:age_at_survey 0.02159 0.02289 1736.27602
site_typeMPA:age_at_survey 0.12217 0.04042 1734.05950
site_typeMPA:kelp_annual_100 0.18151 0.04876 1743.91668
age_at_survey:kelp_annual_100 -0.02467 0.03663 1744.05339
hard_bottom_500:site_typeMPA:age_at_survey 0.34609 0.04690 1739.21871
site_typeMPA:age_at_survey:kelp_annual_100 0.01675 0.04669 1730.99960
t value Pr(>|t|)
(Intercept) 0.720 0.523638
hard_bottom_500 -1.683 0.092580 .
site_typeMPA 6.868 9.01e-12 ***
age_at_survey 0.883 0.377473
kelp_annual_100 -1.528 0.126611
hard_bottom_500:site_typeMPA 9.542 < 2e-16 ***
hard_bottom_500:age_at_survey 0.943 0.345589
site_typeMPA:age_at_survey 3.022 0.002545 **
site_typeMPA:kelp_annual_100 3.722 0.000204 ***
age_at_survey:kelp_annual_100 -0.674 0.500710
hard_bottom_500:site_typeMPA:age_at_survey 7.379 2.45e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100 0.359 0.719863
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500 -0.022
site_typMPA -0.081 -0.016
age_at_srvy 0.016 0.065 -0.055
klp_nnl_100 0.021 0.122 -0.151 0.345
hr__500:_MPA -0.012 -0.304 0.012 0.040 0.014
hrd__500:__ 0.007 -0.006 -0.022 0.123 -0.017 0.008
st_tyMPA:__ -0.004 -0.077 0.025 -0.694 -0.237 -0.055 -0.085
s_MPA:__100 -0.016 -0.125 0.114 -0.283 -0.768 -0.062 0.014
ag_t_:__100 0.033 0.028 -0.220 0.276 0.577 0.027 0.065
h__500:_MPA: -0.006 0.011 -0.037 -0.067 0.011 0.010 -0.481
s_MPA:__:__ -0.024 -0.036 0.239 -0.213 -0.442 -0.064 -0.050
st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500
site_typMPA
age_at_srvy
klp_nnl_100
hr__500:_MPA
hrd__500:__
st_tyMPA:__
s_MPA:__100 0.269
ag_t_:__100 -0.195 -0.451
h__500:_MPA: -0.038 -0.085 -0.033
s_MPA:__:__ 0.161 0.527 -0.761 -0.113 Log+1 Transformed DV
log(B+1) ~ H*ST + K*ST + ST*A (No Scaling)
- Removed year as it had zero variance (singular fit)
- Significant:
H*ST, K*ST, ST*A, ST - Almost significant:
ST - Same conclusions as version without log-transformation (e.g. vs. glm1) except ST is almost significant now
glm5 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescalingattr(glm5, "name") <- "GLMM: log(B+1) ~ H*ST + K*ST + ST*A (Not Scaled)"
summary(glm5)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 *
site_type + site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa)
Data: data_sp
AIC BIC logLik deviance df.resid
-15575.4 -15515.2 7798.7 -15597.4 1737
Scaled residuals:
Min 1Q Median 3Q Max
-2.6738 -0.3367 -0.0312 0.1554 14.2151
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 9.194e-07 0.0009588
bioregion (Intercept) 1.949e-06 0.0013959
Residual 7.571e-06 0.0027516
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.966e-03 9.007e-04 3.624e+00 2.182 0.101554
hard_bottom_500 -9.472e-10 5.746e-10 1.632e+03 -1.648 0.099466
site_typeMPA -3.112e-03 4.109e-04 1.721e+03 -7.575 5.82e-14
kelp_annual_100 -3.866e-08 2.891e-08 1.745e+03 -1.337 0.181362
age_at_survey 1.915e-05 1.765e-05 1.743e+03 1.085 0.278031
hard_bottom_500:site_typeMPA 1.000e-08 1.068e-09 1.574e+03 9.362 < 2e-16
site_typeMPA:kelp_annual_100 1.516e-07 3.598e-08 1.745e+03 4.213 2.64e-05
site_typeMPA:age_at_survey 9.633e-05 2.538e-05 1.734e+03 3.795 0.000153
(Intercept)
hard_bottom_500 .
site_typeMPA ***
kelp_annual_100
age_at_survey
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.232
site_typMPA -0.145 0.297
klp_nnl_100 -0.124 0.130 0.197
age_at_srvy -0.184 0.061 0.361 0.246
h__500:_MPA 0.041 -0.308 -0.740 -0.003 0.033
s_MPA:__100 0.099 -0.124 -0.255 -0.732 -0.200 -0.037
st_tyMPA:__ 0.138 -0.074 -0.521 -0.161 -0.681 -0.050 0.215
fit warnings:
Some predictor variables are on very different scales: consider rescalinglog(B+1) ~ H*ST + K*ST + ST*A (Scaled)
- Year has zero variance (singular fit)
- Significant:
H*ST, K*ST, ST*A - Almost significant: ST
- Same conclusions as version without log-transformation (e.g. vs. glm2) and same conclusions as non-scaled (e.g. vs. glm5)
glm6 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)
attr(glm6, "name") <- "GLMM: log(B+1) ~ H*ST + K*ST + ST*A (Scaled)"
summary(glm6)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type + kelp_annual_100 *
site_type + site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa)
Data: data_sp_scaled
AIC BIC logLik deviance df.resid
4301.3 4361.5 -2139.7 4279.3 1737
Scaled residuals:
Min 1Q Median 3Q Max
-2.6738 -0.3367 -0.0312 0.1554 14.2151
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 0.07978 0.2825
bioregion (Intercept) 0.16910 0.4112
Residual 0.65702 0.8106
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.19720 0.25377 3.04449 0.777 0.493008
hard_bottom_500 -0.04455 0.02703 1632.16947 -1.648 0.099466
site_typeMPA 0.28889 0.04074 1742.21461 7.092 1.92e-12
kelp_annual_100 -0.04526 0.03385 1744.66293 -1.337 0.181362
age_at_survey 0.03038 0.02800 1743.04739 1.085 0.278031
hard_bottom_500:site_typeMPA 0.47034 0.05024 1574.09919 9.362 < 2e-16
site_typeMPA:kelp_annual_100 0.17751 0.04213 1745.08538 4.213 2.64e-05
site_typeMPA:age_at_survey 0.15279 0.04026 1734.05702 3.795 0.000153
(Intercept)
hard_bottom_500 .
site_typeMPA ***
kelp_annual_100
age_at_survey
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 ***
site_typeMPA:age_at_survey ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.024
site_typMPA -0.080 -0.007
klp_nnl_100 0.003 0.130 -0.030
age_at_srvy 0.007 0.061 0.008 0.246
h__500:_MPA -0.013 -0.308 0.026 -0.003 0.033
s_MPA:__100 -0.003 -0.124 -0.021 -0.732 -0.200 -0.037
st_tyMPA:__ 0.003 -0.074 -0.024 -0.161 -0.681 -0.050 0.215 log(B+1) ~ H*ST*A + K*ST*A (Not Scaled)
- Removed year as it had zero variance (singular fit)
- Significant:
H*ST*A, K*ST, ST*A - Almost significant:
H, ST - Same conclusions as version without log-transformation (e.g. vs. glm3)
glm7 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp, REML = FALSE)Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescalingattr(glm7, "name") <- "GLMM: log(B+1) ~ H*ST*A + K*ST*A (Not Scaled)"
summary(glm7)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey +
kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa)
Data: data_sp
AIC BIC logLik deviance df.resid
-15651.1 -15569.1 7840.5 -15681.1 1733
Scaled residuals:
Min 1Q Median 3Q Max
-3.0105 -0.3230 -0.0405 0.1650 14.6102
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 9.051e-07 0.0009514
bioregion (Intercept) 2.093e-06 0.0014466
Residual 7.214e-06 0.0026859
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df
(Intercept) 2.144e-03 9.514e-04 3.941e+00
hard_bottom_500 -1.661e-09 9.417e-10 1.744e+03
site_typeMPA 4.850e-04 6.156e-04 1.739e+03
age_at_survey -2.221e-06 3.145e-05 1.740e+03
kelp_annual_100 -2.039e-08 4.008e-08 1.746e+03
hard_bottom_500:site_typeMPA -1.408e-09 1.854e-09 1.717e+03
hard_bottom_500:age_at_survey 8.557e-11 9.021e-11 1.736e+03
site_typeMPA:age_at_survey -3.384e-04 6.114e-05 1.737e+03
site_typeMPA:kelp_annual_100 1.332e-07 5.319e-08 1.737e+03
age_at_survey:kelp_annual_100 -3.939e-09 5.800e-09 1.744e+03
hard_bottom_500:site_typeMPA:age_at_survey 1.368e-09 1.849e-10 1.739e+03
site_typeMPA:age_at_survey:kelp_annual_100 2.633e-09 7.394e-09 1.731e+03
t value Pr(>|t|)
(Intercept) 2.254 0.0883 .
hard_bottom_500 -1.764 0.0778 .
site_typeMPA 0.788 0.4308
age_at_survey -0.071 0.9437
kelp_annual_100 -0.509 0.6110
hard_bottom_500:site_typeMPA -0.759 0.4479
hard_bottom_500:age_at_survey 0.949 0.3430
site_typeMPA:age_at_survey -5.535 3.58e-08 ***
site_typeMPA:kelp_annual_100 2.504 0.0124 *
age_at_survey:kelp_annual_100 -0.679 0.4972
hard_bottom_500:site_typeMPA:age_at_survey 7.400 2.10e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100 0.356 0.7218
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500 -0.315
site_typMPA -0.174 0.398
age_at_srvy -0.290 0.669 0.428
klp_nnl_100 -0.130 0.118 0.185 0.309
hr__500:_MPA 0.118 -0.431 -0.889 -0.318 -0.053
hrd__500:__ 0.233 -0.803 -0.353 -0.813 -0.094 0.405
st_tyMPA:__ 0.154 -0.350 -0.830 -0.501 -0.155 0.736 0.413
s_MPA:__100 0.093 -0.088 -0.179 -0.229 -0.714 -0.040 0.068
ag_t_:__100 0.054 -0.035 -0.108 -0.248 -0.708 0.043 0.065
h__500:_MPA: -0.115 0.391 0.740 0.387 0.050 -0.826 -0.481
s_MPA:__:__ -0.037 0.018 0.069 0.187 0.536 0.058 -0.050
st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500
site_typMPA
age_at_srvy
klp_nnl_100
hr__500:_MPA
hrd__500:__
st_tyMPA:__
s_MPA:__100 0.122
ag_t_:__100 0.123 0.530
h__500:_MPA: -0.897 0.065 -0.033
s_MPA:__:__ -0.060 -0.748 -0.761 -0.113
fit warnings:
Some predictor variables are on very different scales: consider rescalinglog(B+1) ~ H*ST*A + K*ST*A (Scaled)
- Removed year as it had zero variance (singular fit)
- Significant:
H*ST*A, H*ST, K*ST, ST*A, ST - Same conclusions as version without log-transformation (e.g. vs. glm4) but different from non-scaled version b/c ST now highly significant on its own, and H*ST is significant on its own, hard bottom on its own is nearly significant
glm8 <- lmer(log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa), data = data_sp_scaled, REML = FALSE)
attr(glm8, "name") <- "GLMM: log(B+1) ~ H*ST*A + K*ST*A (Scaled)"
summary(glm8)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2 ~ hard_bottom_500 * site_type * age_at_survey +
kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa)
Data: data_sp_scaled
AIC BIC logLik deviance df.resid
4225.6 4307.6 -2097.8 4195.6 1733
Scaled residuals:
Min 1Q Median 3Q Max
-3.0105 -0.3230 -0.0405 0.1650 14.6102
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 0.07854 0.2803
bioregion (Intercept) 0.18158 0.4261
Residual 0.62600 0.7912
Number of obs: 1748, groups: affiliated_mpa, 23; bioregion, 3
Fixed effects:
Estimate Std. Error df
(Intercept) 0.18905 0.26173 3.00929
hard_bottom_500 -0.04457 0.02642 1637.48328
site_typeMPA 0.28265 0.04105 1741.42207
age_at_survey 0.02536 0.02861 1743.62208
kelp_annual_100 -0.06233 0.04056 1743.89705
hard_bottom_500:site_typeMPA 0.47054 0.04920 1587.09477
hard_bottom_500:age_at_survey 0.02167 0.02285 1736.20266
site_typeMPA:age_at_survey 0.12248 0.04035 1734.00634
site_typeMPA:kelp_annual_100 0.18165 0.04868 1743.85438
age_at_survey:kelp_annual_100 -0.02483 0.03656 1744.00838
hard_bottom_500:site_typeMPA:age_at_survey 0.34648 0.04682 1739.16854
site_typeMPA:age_at_survey:kelp_annual_100 0.01660 0.04661 1730.96135
t value Pr(>|t|)
(Intercept) 0.722 0.522162
hard_bottom_500 -1.687 0.091773 .
site_typeMPA 6.885 8.02e-12 ***
age_at_survey 0.886 0.375603
kelp_annual_100 -1.537 0.124554
hard_bottom_500:site_typeMPA 9.564 < 2e-16 ***
hard_bottom_500:age_at_survey 0.949 0.342965
site_typeMPA:age_at_survey 3.035 0.002439 **
site_typeMPA:kelp_annual_100 3.732 0.000196 ***
age_at_survey:kelp_annual_100 -0.679 0.497184
hard_bottom_500:site_typeMPA:age_at_survey 7.400 2.10e-13 ***
site_typeMPA:age_at_survey:kelp_annual_100 0.356 0.721838
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500 -0.022
site_typMPA -0.080 -0.016
age_at_srvy 0.016 0.065 -0.055
klp_nnl_100 0.021 0.122 -0.151 0.345
hr__500:_MPA -0.012 -0.304 0.012 0.040 0.014
hrd__500:__ 0.007 -0.006 -0.022 0.123 -0.017 0.008
st_tyMPA:__ -0.004 -0.077 0.025 -0.694 -0.237 -0.055 -0.085
s_MPA:__100 -0.016 -0.125 0.114 -0.283 -0.768 -0.063 0.014
ag_t_:__100 0.033 0.028 -0.220 0.276 0.577 0.027 0.065
h__500:_MPA: -0.006 0.010 -0.037 -0.067 0.011 0.010 -0.481
s_MPA:__:__ -0.024 -0.036 0.239 -0.213 -0.442 -0.064 -0.050
st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500
site_typMPA
age_at_srvy
klp_nnl_100
hr__500:_MPA
hrd__500:__
st_tyMPA:__
s_MPA:__100 0.269
ag_t_:__100 -0.195 -0.451
h__500:_MPA: -0.038 -0.085 -0.033
s_MPA:__:__ 0.161 0.527 -0.761 -0.113 Log + Small Constant Transformation
log(B+c) ~ H*ST + K*ST + ST*A (No Scaling)
- Now there is nonzero variance for year
H*ST and K*STis still significant butST*Ais not anymore- Now intercept is significant, H alone is significant, Age alone is significant
- Changes compared to the non-transformed and log+1 transformation!
glm9 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp, REML = FALSE)Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescalingattr(glm9, "name") <- "GLMM: log(B+c) ~ H*ST + K*ST + ST*A (Not Scaled)"
summary(glm9)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 *
site_type + site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa) + (1 | year)
Data: data_sp
AIC BIC logLik deviance df.resid
9139.7 9205.3 -4557.9 9115.7 1736
Scaled residuals:
Min 1Q Median 3Q Max
-3.2706 -0.5521 -0.0645 0.5798 3.2244
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 2.3174 1.5223
year (Intercept) 0.1354 0.3679
bioregion (Intercept) 8.0802 2.8426
Residual 10.2866 3.2073
Number of obs: 1748, groups: affiliated_mpa, 23; year, 21; bioregion, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.213e+01 1.737e+00 3.365e+00 -6.986 0.004100
hard_bottom_500 -1.568e-06 6.769e-07 1.698e+03 -2.316 0.020661
site_typeMPA -1.091e+00 4.817e-01 1.723e+03 -2.265 0.023641
kelp_annual_100 -4.279e-05 3.403e-05 1.746e+03 -1.257 0.208803
age_at_survey 1.008e-01 2.487e-02 5.776e+01 4.054 0.000152
hard_bottom_500:site_typeMPA 5.806e-06 1.259e-06 1.667e+03 4.610 4.33e-06
site_typeMPA:kelp_annual_100 1.386e-04 4.210e-05 1.730e+03 3.291 0.001017
site_typeMPA:age_at_survey -2.798e-02 2.963e-02 1.716e+03 -0.944 0.345069
(Intercept) **
hard_bottom_500 *
site_typeMPA *
kelp_annual_100
age_at_survey ***
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 **
site_typeMPA:age_at_survey
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.142
site_typMPA -0.086 0.285
klp_nnl_100 -0.075 0.127 0.192
age_at_srvy -0.111 0.048 0.298 0.203
h__500:_MPA 0.022 -0.291 -0.742 0.002 0.029
s_MPA:__100 0.060 -0.126 -0.251 -0.728 -0.166 -0.042
st_tyMPA:__ 0.084 -0.075 -0.516 -0.161 -0.566 -0.052 0.216
fit warnings:
Some predictor variables are on very different scales: consider rescalinglog(B+c) ~ H*ST + K*ST + ST*A (Scaled)
- Now there is nonzero variance for year
H*ST and K*STis still significant butST*Ais not anymore- Now H alone is significant, Age alone is significant, ST is significant
- Changes compared to the non-transformed and log+1 transformation, similar to the non-scaled log+c version (e.g. glmm #9) except intercept is not significant in this version
glm10 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp_scaled, REML = FALSE)
attr(glm10, "name") <- "GLMM: log(B+c) ~ H*ST + K*ST + ST*A (Scaled)"
summary(glm10)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 *
site_type + site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa) + (1 | year)
Data: data_sp_scaled
AIC BIC logLik deviance df.resid
3935.5 4001.1 -1955.7 3911.5 1736
Scaled residuals:
Min 1Q Median 3Q Max
-3.2706 -0.5521 -0.0645 0.5798 3.2244
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 0.118046 0.34358
year (Intercept) 0.006895 0.08304
bioregion (Intercept) 0.411906 0.64180
Residual 0.523924 0.72383
Number of obs: 1748, groups: affiliated_mpa, 23; year, 21; bioregion, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.39905 0.38581 3.15169 1.034 0.373731
hard_bottom_500 -0.05650 0.02439 1698.09694 -2.316 0.020659
site_typeMPA 0.15380 0.03644 1718.90196 4.221 2.56e-05
kelp_annual_100 -0.03838 0.03053 1745.56360 -1.257 0.208789
age_at_survey 0.12251 0.03022 57.75493 4.054 0.000153
hard_bottom_500:site_typeMPA 0.20922 0.04538 1666.64841 4.610 4.33e-06
site_typeMPA:kelp_annual_100 0.12428 0.03776 1729.73277 3.291 0.001017
site_typeMPA:age_at_survey -0.03400 0.03600 1716.16888 -0.944 0.345069
(Intercept)
hard_bottom_500 *
site_typeMPA ***
kelp_annual_100
age_at_survey ***
hard_bottom_500:site_typeMPA ***
site_typeMPA:kelp_annual_100 **
site_typeMPA:age_at_survey
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA k__100 ag_t_s h__500: s_MPA:__1
hrd_btt_500 -0.015
site_typMPA -0.047 -0.007
klp_nnl_100 0.001 0.127 -0.030
age_at_srvy 0.021 0.048 0.007 0.203
h__500:_MPA -0.008 -0.291 0.026 0.002 0.029
s_MPA:__100 -0.001 -0.126 -0.022 -0.728 -0.166 -0.042
st_tyMPA:__ 0.002 -0.075 -0.025 -0.161 -0.566 -0.052 0.216 log(B+c) ~ H*ST*A + K*ST*A (Not Scaled)
- Now there is nonzero variance for year
H*ST*Ais still significant but less so;ST*A and K*STis significant- Intercept, H alone are now both highly significant; ST no longer significant
- Now significant:
H*A - Changes compared to the non-transformed and log+1 transformation.
glm11 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp, REML = FALSE)Warning: Some predictor variables are on very different scales: consider
rescaling
Warning: Some predictor variables are on very different scales: consider
rescalingattr(glm11, "name") <- "GLMM: log(B+c) ~ H*ST*A + K*ST*A (Not Scaled)"
summary(glm11)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey +
kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa) + (1 | year)
Data: data_sp
AIC BIC logLik deviance df.resid
9129.2 9216.7 -4548.6 9097.2 1732
Scaled residuals:
Min 1Q Median 3Q Max
-3.1972 -0.5341 -0.0619 0.5650 3.3157
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 2.3012 1.5170
year (Intercept) 0.1151 0.3392
bioregion (Intercept) 8.1154 2.8487
Residual 10.1863 3.1916
Number of obs: 1748, groups: affiliated_mpa, 23; year, 21; bioregion, 3
Fixed effects:
Estimate Std. Error df
(Intercept) -1.162e+01 1.759e+00 3.509e+00
hard_bottom_500 -3.469e-06 1.129e-06 1.746e+03
site_typeMPA 1.612e-01 7.352e-01 1.733e+03
age_at_survey 4.065e-02 3.939e-02 3.986e+02
kelp_annual_100 -3.321e-05 4.814e-05 1.733e+03
hard_bottom_500:site_typeMPA 1.732e-06 2.220e-06 1.732e+03
hard_bottom_500:age_at_survey 2.276e-07 1.080e-07 1.727e+03
site_typeMPA:age_at_survey -1.825e-01 7.280e-02 1.722e+03
site_typeMPA:kelp_annual_100 1.637e-04 6.335e-05 1.724e+03
age_at_survey:kelp_annual_100 -2.962e-06 6.952e-06 1.739e+03
hard_bottom_500:site_typeMPA:age_at_survey 4.974e-07 2.204e-07 1.730e+03
site_typeMPA:age_at_survey:kelp_annual_100 -3.620e-06 8.800e-06 1.718e+03
t value Pr(>|t|)
(Intercept) -6.606 0.00427 **
hard_bottom_500 -3.072 0.00216 **
site_typeMPA 0.219 0.82644
age_at_survey 1.032 0.30261
kelp_annual_100 -0.690 0.49031
hard_bottom_500:site_typeMPA 0.780 0.43548
hard_bottom_500:age_at_survey 2.109 0.03513 *
site_typeMPA:age_at_survey -2.507 0.01227 *
site_typeMPA:kelp_annual_100 2.583 0.00987 **
age_at_survey:kelp_annual_100 -0.426 0.67009
hard_bottom_500:site_typeMPA:age_at_survey 2.257 0.02413 *
site_typeMPA:age_at_survey:kelp_annual_100 -0.411 0.68083
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500 -0.203
site_typMPA -0.111 0.391
age_at_srvy -0.186 0.628 0.403
klp_nnl_100 -0.084 0.115 0.181 0.296
hr__500:_MPA 0.074 -0.422 -0.890 -0.298 -0.050
hrd__500:__ 0.149 -0.802 -0.351 -0.765 -0.093 0.402
st_tyMPA:__ 0.099 -0.347 -0.829 -0.473 -0.152 0.735 0.411
s_MPA:__100 0.060 -0.087 -0.175 -0.219 -0.711 -0.044 0.067
ag_t_:__100 0.035 -0.032 -0.107 -0.239 -0.709 0.042 0.063
h__500:_MPA: -0.074 0.388 0.740 0.364 0.047 -0.825 -0.479
s_MPA:__:__ -0.024 0.016 0.067 0.179 0.533 0.060 -0.048
st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500
site_typMPA
age_at_srvy
klp_nnl_100
hr__500:_MPA
hrd__500:__
st_tyMPA:__
s_MPA:__100 0.120
ag_t_:__100 0.122 0.527
h__500:_MPA: -0.897 0.068 -0.032
s_MPA:__:__ -0.058 -0.747 -0.757 -0.115
fit warnings:
Some predictor variables are on very different scales: consider rescalinglog(B+c) ~ H*ST*A + K*ST*A (Scaled)
- Now there is nonzero variance for year
- Significant:
H*ST*A (but less so), H*ST, K*ST; now: H, ST, A, H*A - Changes compared to the non-transformed and log+1 transformation.
ST*Ano longer significant, andH, ST, A, H*Aare now significant on their own.
glm12 <- lmer(log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey + kelp_annual_100 * site_type * age_at_survey +
(1|bioregion) + (1|affiliated_mpa) + (1|year), data = data_sp_scaled, REML = FALSE)
attr(glm12, "name") <- "GLMM: log(B+c) ~ H*ST*A + K*ST*A (Scaled)"
summary(glm12)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: log_kg_per_m2_c ~ hard_bottom_500 * site_type * age_at_survey +
kelp_annual_100 * site_type * age_at_survey + (1 | bioregion) +
(1 | affiliated_mpa) + (1 | year)
Data: data_sp_scaled
AIC BIC logLik deviance df.resid
3925.0 4012.4 -1946.5 3893.0 1732
Scaled residuals:
Min 1Q Median 3Q Max
-3.1972 -0.5341 -0.0619 0.5650 3.3157
Random effects:
Groups Name Variance Std.Dev.
affiliated_mpa (Intercept) 0.117206 0.34235
year (Intercept) 0.005862 0.07656
bioregion (Intercept) 0.413363 0.64293
Residual 0.518815 0.72029
Number of obs: 1748, groups: affiliated_mpa, 23; year, 21; bioregion, 3
Fixed effects:
Estimate Std. Error df
(Intercept) 0.39457 0.38632 3.14695
hard_bottom_500 -0.05657 0.02429 1698.00635
site_typeMPA 0.14844 0.03743 1717.60963
age_at_survey 0.12560 0.03047 66.06507
kelp_annual_100 -0.05195 0.03720 1743.96625
hard_bottom_500:site_typeMPA 0.21192 0.04527 1669.08203
hard_bottom_500:age_at_survey 0.04417 0.02095 1726.98524
site_typeMPA:age_at_survey -0.04867 0.03677 1715.10015
site_typeMPA:kelp_annual_100 0.11971 0.04443 1726.04778
age_at_survey:kelp_annual_100 -0.01431 0.03358 1739.06377
hard_bottom_500:site_typeMPA:age_at_survey 0.09651 0.04276 1729.97794
site_typeMPA:age_at_survey:kelp_annual_100 -0.01748 0.04250 1718.08058
t value Pr(>|t|)
(Intercept) 1.021 0.379075
hard_bottom_500 -2.329 0.019983 *
site_typeMPA 3.966 7.60e-05 ***
age_at_survey 4.122 0.000107 ***
kelp_annual_100 -1.397 0.162692
hard_bottom_500:site_typeMPA 4.681 3.08e-06 ***
hard_bottom_500:age_at_survey 2.109 0.035127 *
site_typeMPA:age_at_survey -1.323 0.185850
site_typeMPA:kelp_annual_100 2.694 0.007123 **
age_at_survey:kelp_annual_100 -0.426 0.670088
hard_bottom_500:site_typeMPA:age_at_survey 2.257 0.024134 *
site_typeMPA:age_at_survey:kelp_annual_100 -0.411 0.680836
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hr__500 st_MPA ag_t_s k__100 hr__500:_MPA h__500:__
hrd_btt_500 -0.014
site_typMPA -0.050 -0.016
age_at_srvy 0.024 0.053 -0.046
klp_nnl_100 0.013 0.121 -0.150 0.292
hr__500:_MPA -0.008 -0.288 0.012 0.036 0.018
hrd__500:__ 0.005 -0.008 -0.021 0.124 -0.020 0.009
st_tyMPA:__ -0.002 -0.078 0.024 -0.597 -0.236 -0.058 -0.086
s_MPA:__100 -0.010 -0.126 0.112 -0.239 -0.764 -0.066 0.016
ag_t_:__100 0.020 0.030 -0.219 0.233 0.575 0.029 0.063
h__500:_MPA: -0.004 0.009 -0.037 -0.065 0.010 0.006 -0.479
s_MPA:__:__ -0.014 -0.038 0.239 -0.179 -0.440 -0.063 -0.048
st_MPA:__ s_MPA:__1 a__:__ h__500:_MPA:
hrd_btt_500
site_typMPA
age_at_srvy
klp_nnl_100
hr__500:_MPA
hrd__500:__
st_tyMPA:__
s_MPA:__100 0.269
ag_t_:__100 -0.194 -0.448
h__500:_MPA: -0.036 -0.084 -0.032
s_MPA:__:__ 0.161 0.526 -0.757 -0.115 Quick overview:
- Not much difference between scaled and non-scaled versions in terms of conclusions
- Not much difference between non-transformed and log+1 transformed versions in terms of conclusions
- Clear differences in the log+c transformed versions compared to the other two.
- Next: Need to see whether scaling or the transformations yields better fit.
Comparison & Diagnostics
model_fit <- bind_rows(model_performance(glm1, estimator = "ML") %>% mutate(Model = attr(glm1, "name")),
model_performance(glm2, estimator = "ML") %>% mutate(Model = attr(glm2, "name")),
model_performance(glm3, estimator = "ML") %>% mutate(Model = attr(glm3, "name")),
model_performance(glm4, estimator = "ML") %>% mutate(Model = attr(glm4, "name"))) %>%
mutate(deltaAIC = AIC - min(AIC, na.rm = T)) %>%
dplyr::select(Model, everything()) %>%
arrange(AIC)
print(model_fit)# Indices of model performance
Model | AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | deltaAIC
----------------------------------------------------------------------------------------------------------------------------------------------
GLMM: B ~ H*ST*A + K*ST*A (Not Scaled) | -15611.275 | -15610.998 | -15529.282 | 0.382 | 0.129 | 0.291 | 0.003 | 0.003 | 0.000
GLMM: B ~ H*ST + K*ST + ST*A (Not Scaled) | -15536.055 | -15535.903 | -15475.927 | 0.345 | 0.099 | 0.272 | 0.003 | 0.003 | 75.220
GLMM: B ~ H*ST*A + K*ST*A (Scaled) | 4231.592 | 4231.869 | 4313.586 | 0.382 | 0.129 | 0.291 | 0.788 | 0.793 | 19842.867
GLMM: B ~ H*ST + K*ST + ST*A (Scaled) | 4306.812 | 4306.964 | 4366.941 | 0.345 | 0.099 | 0.272 | 0.807 | 0.812 | 19918.087model_fit2 <- bind_rows(model_performance(glm5, estimator = "ML") %>% mutate(Model = attr(glm5, "name")),
model_performance(glm6, estimator = "ML") %>% mutate(Model = attr(glm6, "name")),
model_performance(glm7, estimator = "ML") %>% mutate(Model = attr(glm7, "name")),
model_performance(glm8, estimator = "ML") %>% mutate(Model = attr(glm8, "name"))) %>%
mutate(deltaAIC = AIC - min(AIC, na.rm = T)) %>%
dplyr::select(Model, everything()) %>%
arrange(AIC)
print(model_fit2)# Indices of model performance
Model | AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | deltaAIC
-----------------------------------------------------------------------------------------------------------------------------------------------------
GLMM: log(B+1) ~ H*ST*A + K*ST*A (Not Scaled) | -15651.052 | -15650.775 | -15569.059 | 0.385 | 0.129 | 0.294 | 0.003 | 0.003 | 0.000
GLMM: log(B+1) ~ H*ST + K*ST + ST*A (Not Scaled) | -15575.372 | -15575.220 | -15515.243 | 0.347 | 0.100 | 0.275 | 0.003 | 0.003 | 75.680
GLMM: log(B+1) ~ H*ST*A + K*ST*A (Scaled) | 4225.648 | 4225.925 | 4307.641 | 0.385 | 0.129 | 0.294 | 0.787 | 0.791 | 19876.700
GLMM: log(B+1) ~ H*ST + K*ST + ST*A (Scaled) | 4301.328 | 4301.480 | 4361.457 | 0.347 | 0.100 | 0.275 | 0.806 | 0.811 | 19952.380model_fit3 <- bind_rows(model_performance(glm9, estimator = "ML") %>% mutate(Model = attr(glm9, "name")),
model_performance(glm10, estimator = "ML") %>% mutate(Model = attr(glm10, "name")),
model_performance(glm11, estimator = "ML") %>% mutate(Model = attr(glm11, "name")),
model_performance(glm12 ,estimator = "ML") %>% mutate(Model = attr(glm12, "name"))) %>%
mutate(deltaAIC = AIC - min(AIC, na.rm = T)) %>%
dplyr::select(Model, everything()) %>%
arrange(AIC)
print(model_fit3)# Indices of model performance
Model | AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | deltaAIC
----------------------------------------------------------------------------------------------------------------------------------------------
GLMM: log(B+c) ~ H*ST*A + K*ST*A (Scaled) | 3924.961 | 3925.276 | 4012.421 | 0.526 | 0.035 | 0.508 | 0.714 | 0.720 | 0.000
GLMM: log(B+c) ~ H*ST + K*ST + ST*A (Scaled) | 3935.475 | 3935.655 | 4001.070 | 0.521 | 0.030 | 0.506 | 0.717 | 0.724 | 10.514
GLMM: log(B+c) ~ H*ST*A + K*ST*A (Not Scaled) | 9129.190 | 9129.505 | 9216.650 | 0.526 | 0.035 | 0.508 | 3.164 | 3.192 | 5204.229
GLMM: log(B+c) ~ H*ST + K*ST + ST*A (Not Scaled) | 9139.704 | 9139.884 | 9205.299 | 0.521 | 0.030 | 0.506 | 3.179 | 3.207 | 5214.743par(mfrow = c(2, 2))
plot(residuals(glm1) ~ fitted(glm1), main = attr(glm1, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm1), residuals(glm1)), col = "blue", lwd = 2)
plot(residuals(glm2) ~ fitted(glm2), main = attr(glm2, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm2), residuals(glm2)), col = "blue", lwd = 2)
plot(residuals(glm3) ~ fitted(glm3), main = attr(glm3, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm3), residuals(glm3)), col = "blue", lwd = 2)
plot(residuals(glm4) ~ fitted(glm4), main = attr(glm4, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm4), residuals(glm4)), col = "blue", lwd = 2)
par(mfrow = c(2, 2))
plot(residuals(glm5) ~ fitted(glm5), main = attr(glm5, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm5), residuals(glm5)), col = "blue", lwd = 2)
plot(residuals(glm6) ~ fitted(glm6), main = attr(glm6, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm6), residuals(glm6)), col = "blue", lwd = 2)
plot(residuals(glm7) ~ fitted(glm7), main = attr(glm7, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm7), residuals(glm7)), col = "blue", lwd = 2)
plot(residuals(glm8) ~ fitted(glm8), main = attr(glm8, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm8), residuals(glm8)), col = "blue", lwd = 2)
par(mfrow = c(2, 2))
plot(residuals(glm9) ~ fitted(glm9), main = attr(glm9, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm9), residuals(glm9)), col = "blue", lwd = 2)
plot(residuals(glm10) ~ fitted(glm10), main = attr(glm10, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm10), residuals(glm10)), col = "blue", lwd = 2)
plot(residuals(glm11) ~ fitted(glm11), main = attr(glm11, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm11), residuals(glm11)), col = "blue", lwd = 2)
plot(residuals(glm12) ~ fitted(glm12), main = attr(glm12, "name"))
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted(glm12), residuals(glm12)), col = "blue", lwd = 2)
par(mfrow = c(2, 2))
qqnorm(residuals(glm1), main = paste(attr(glm1, "name"), "QQ Plot")); qqline(residuals(glm1))
qqnorm(residuals(glm2), main = paste(attr(glm2, "name"), "QQ Plot")); qqline(residuals(glm2))
qqnorm(residuals(glm3), main = paste(attr(glm3, "name"), "QQ Plot")); qqline(residuals(glm3))
qqnorm(residuals(glm4), main = paste(attr(glm4, "name"), "QQ Plot")); qqline(residuals(glm4))
par(mfrow = c(2, 2))
qqnorm(residuals(glm5), main = paste(attr(glm5, "name"), "QQ Plot")); qqline(residuals(glm5))
qqnorm(residuals(glm6), main = paste(attr(glm6, "name"), "QQ Plot")); qqline(residuals(glm6))
qqnorm(residuals(glm7), main = paste(attr(glm7, "name"), "QQ Plot")); qqline(residuals(glm7))
qqnorm(residuals(glm8), main = paste(attr(glm8, "name"), "QQ Plot")); qqline(residuals(glm8))
par(mfrow = c(2, 2))
qqnorm(residuals(glm9), main = paste(attr(glm9, "name"), "QQ Plot")); qqline(residuals(glm9))
qqnorm(residuals(glm10), main = paste(attr(glm10, "name"), "QQ Plot")); qqline(residuals(glm10))
qqnorm(residuals(glm11), main = paste(attr(glm11, "name"), "QQ Plot")); qqline(residuals(glm11))
qqnorm(residuals(glm12), main = paste(attr(glm12, "name"), "QQ Plot")); qqline(residuals(glm12))
Partial Residual Plots
B ~ H*ST + K*ST + ST*A (Not Scaled)
plot(predictorEffects(glm1, ~ age_at_survey, partial.residuals = TRUE),
id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm1, ~ hard_bottom_500, partial.residuals = TRUE),
id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm1, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)),
id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)
# Error likely means that the extremes in the hard bottom value are causing challenges?
# Looks like large number of low values in kelp canopy cover .... interesting.B ~ H*ST + K*ST + ST*A (Scaled)
plot(predictorEffects(glm2, ~ age_at_survey, partial.residuals = TRUE),
id = list(n = 1), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm2, ~ hard_bottom_500, partial.residuals = TRUE),
id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm2, ~ kelp_annual_100, partial.residuals = TRUE),
id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
B ~ H*ST*A + K*ST*A (Not Scaled)
plot(predictorEffects(glm3, ~ age_at_survey, partial.residuals = TRUE),
axes = list(x = list(cex = 0.8), y = list(cex = 0.8, lim = c(-0.001, 0.02))),
lattice = list(strip = list(cex = 0.5)),
id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm3, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm3, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 35)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
B ~ H*ST*A + K*ST*A (Scaled)
plot(predictorEffects(glm4, ~ age_at_survey, partial.residuals = TRUE),
axes = list( x = list(rotate = 25, cex = 0.8), y = list(cex = 0.8)),
lattice = list(strip = list(cex = 0.5)),
id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm4, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25), y = list(lim = c(-1, 5))), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm4, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25), y = list(lim = c(-1, 5))), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
log(B+1) ~ H*ST + K*ST + ST*A (Not Scaled)
plot(predictorEffects(glm5, ~ age_at_survey, partial.residuals = TRUE),
id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm5, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm5, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
# Error likely means that the extremes in the hard bottom value are causing challenges?log(B+1) ~ H*ST + K*ST + ST*A (Scaled)
plot(predictorEffects(glm6, ~ age_at_survey, partial.residuals = TRUE),
id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm6, ~ hard_bottom_500, partial.residuals = TRUE), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm6, ~ kelp_annual_100, partial.residuals = TRUE), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
log(B+1) ~ H*ST*A + K*ST*A (Not Scaled)
plot(predictorEffects(glm7, ~ age_at_survey, partial.residuals = TRUE),
axes = list( x = list(cex = 0.8), y = list(cex = 0.8)),
lattice = list(strip = list(cex = 0.5)),
id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm7, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm7, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
log(B+1) ~ H*ST*A + K*ST*A (Scaled)
plot(predictorEffects(glm8, ~ age_at_survey, partial.residuals = TRUE),
axes = list( x = list(rotate = 25, cex = 0.8), y = list(cex = 0.8)),
lattice = list(strip = list(cex = 0.5)),
id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm8, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm8, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
log(B+c) ~ H*ST + K*ST + ST*A (Not Scaled)
plot(predictorEffects(glm9, ~ age_at_survey, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm9, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm9, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
log(B+c) ~ H*ST + K*ST + ST*A (Scaled)
plot(predictorEffects(glm10, ~ age_at_survey, partial.residuals = TRUE),
id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm10, ~ hard_bottom_500, partial.residuals = TRUE),
id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm10, ~ kelp_annual_100, partial.residuals = TRUE), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
log(B+c) ~ H*ST*A + K*ST*A (Not Scaled)
plot(predictorEffects(glm11, ~ age_at_survey, partial.residuals = TRUE),
axes = list( x = list(cex = 0.8), y = list(cex = 0.8)),
lattice = list(strip = list(cex = 0.5)),
id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm11, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm11, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25)), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
log(B+c) ~ H*ST*A + K*ST*A (Scaled)
plot(predictorEffects(glm12, ~ age_at_survey, partial.residuals = TRUE),
axes = list( x = list(cex = 0.8), y = list(cex = 0.8)),
lattice = list(strip = list(cex = 0.5)),
id = list(n = 1, cex = 0.5), residuals.pch = 19, residuals.cex = 0.2)
plot(predictorEffects(glm12, ~ hard_bottom_500, partial.residuals = TRUE),
axes = list(x = list(rotate = 25), y = list(lim = c(-1, 2))), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
plot(predictorEffects(glm12, ~ kelp_annual_100, partial.residuals = TRUE),
axes = list(x = list(rotate = 25), y = list(lim = c(-1, 2))), id = list(n = 1),
residuals.pch = 19,
residuals.cex = 0.2)
GLMM Summary
- The diagnostic plots don’t show drastic differences among the non-transformed vs. log+1 transformed. Within those two, the non-scaled, 3-way interaction tends to have the lowest AIC.
- The log+c transformation causes challenges because the zeroes become negative values (log of a value less than 1 is negative).
- Scaling the kelp canopy cover may be creating challenges (scaled values go from -0.5 to 6… suggests lots of zeroes and some very high values). This could be because kelp is scaled across all years, rather than within each year?
Gamma GLMs
1. Fit Gamma GLM + 2way
# gglm1 <- glmer(kg_per_m2_c ~ hard_bottom_500 * site_type + kelp_annual_100 * site_type + site_type * age_at_survey +
# (1|bioregion) + (1|affiliated_mpa) + (1|year),
# family = Gamma(link = "log"), data = data_sp)
# does not fit